heterozygosity <- function(inbfile=NULL,hist_breaks=NULL,outputfile="heterozygosity"){

	diff.het <- function(data){
		res <- as.numeric(data[c(4)]) - as.numeric(data[c(3)])
	}

	if(is.null(hist_breaks)){
		rm(hist_breaks)	
	}

	# creation d'un nombre aléatoire servant de prefixe
	randnb <- round(runif(10,min=0,max=1000))
	prefix <- randnb[1]

	system(paste("plink --noweb --bfile ",inbfile," --het --out ",prefix,"het",sep=""))
	het <- read.table(file=paste(prefix,"het.het",sep=""),h=T)
	het$diffEO <- apply(X=het,MARGIN=1,FUN=diff.het)
	
	write.table(het,file=paste("summary-het",outputfile,".txt",sep=""),quote=F,row.names=F)	

	png(paste("plots-",outputfile,".png",sep=""),height=580,width=920)
	par(mfrow=c(1,2))
	plot(het$diffEO)
	hist(het$diffEO,breaks=hist_breaks)
	dev.off()

	cat(paste("\n see plots-",outputfile,".png and summary-het",outputfile,".txt \n",sep=""))
}
